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The  paper  focuses  on  numerical  simulation  of  neutral  entrainment  in  a  Field  Reversed 
Configuration  thruster  with  a  neon  propellant.  A  fully  kinetic  approach  based  on  a  Particle- 
In-Cell  method  of  modeling  plasmas  and  the  Direct  Simulation  Monte  Carlo  method  of 
modeling  neutrals  is  applied.  An  implicit  PIC  code  Celeste3D  extended  to  compute  neu¬ 
tral  transport  and  collision  processes  between  neutral  and  charged  particles  is  used.  The 
current  collision  model  includes  all  key  elastic,  collisional  radiative,  and  reaction  processes 
in  neon.  Temporal  evolution  of  gas  and  plasma  macroparameters  and  thruster  performance 
properties  is  presented. 


I.  Introduction 

High-power  electric  propulsion  is  an  enabling  technology  for  future  space  missions,  but  one  that  also 
presents  a  number  of  technical  challenges.  First,  it  is  critically  important  to  operate  a  thruster  at  very 
high  efficiency,  otherwise  the  system  limitations  due  to  heat  rejection  become  insurmountable.  Second,  high 
power  can  be  efficiently  delivered  into  a  plasma  by  raising  its  temperature  (provided  radiative  losses  are  not 
excessive),  i.e.  its  specific  energy;  consequently,  the  specific  impulse  can  take  large  values,  even  in  the  excess 
of  10,000  sec.  The  optimal  conditions  from  the  standpoint  of  low  radiation  loss,  acceptable  ionization  cost, 
and  confinement  requirements,  are  usually  achieved  for  low-Z  plasma  at  high  (50-100  eV)  temperature  and 
in  strong  magnetic  fields.  These  conditions  are  typically  obtained  in  Field  Reversed  Configuration  (FRC) 
(see,  for  example,  Ref.  1)  plasma.  The  FRC  is  a  self-organized  magnetized  plasma  structure  in  the  shape 
of  a  highly  compact  toroid.  The  magnetic  field  is  mostly  in  the  poloidal  direction,2  generated  by  internal 
(toroidal)  currents.  The  ratio  of  plasma  pressure  to  magnetic  pressure  is  close  to  unity,  i.e.  the  highest 
plasma  density  that  can  be  attained  for  given  external  magnets.  The  poloidal  field  also  contributes  to  the 
particle  confinement.  Starting  from  a  background  uniform  plasma  at  a  constant  axial  (bias)  field,  the  FRC 
can  be  formed  by  pulsing  external  coils  and  reversing  the  applied  field,  inducing  currents  at  the  plasma 
boundary  and  "pinching”  the  plasma  at  both  ends  (see  Fig.  1).  The  initial  bias  field  is  trapped  inside  the 
plasma  and  forced  to  reconnect  at  the  end  points  (separatrix),  creating  an  elongated  toroidal  shape. 

In  contrast  to  some  other  high  Isp  plasma  propulsion  concepts,3  the  FRC  is  completely  magnetically 
insulated  from  the  external  field:  the  plasma  is  not  tied  to  an  external  field  line,  and  the  FRC  can  readily 
detach  from  the  confining  external  field.  It  can  also  be  translated  and  accelerated  by  applying  a  gradient  of 
magnetic  pressure  using  pulsed  external  coils.  The  FRC  can  therefore  be  efficiently  accelerated  to  provide 
thrust,  operates  at  a  temperature  that  is  optimal  for  ionization  and  is  well  confined.  The  basic  concept  of 
operations  has  been  demonstrated4  and  more  recent  research  has  led  to  further  optimization  of  the  formation 
process,5  while  methods  for  increased  efficiency  through  energy  recovery  in  the  electrical  circuits  are  currently 
being  developed.  One  of  the  latest  design  iterations6  provides  plasma  velocities  in  the  10-40  km/sec  range, 
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Figure  1.  Schematic  of  basic  FRC  formation  by  field  reversal. 


a  desirable  regime  for  USAF  applications.  However,  many  aspects  of  the  flow  development  and  device 
operation,  such  as  power  and  mass  utilization  efficiencies,  still  remain  to  be  determined  more  precisely. 

One  of  the  key  features  of  a  versatile  FRC  thruster  is  its  ability  to  increase  thrust  level  and  operational 
efficiency  through  the  process  of  entrainment  of  neutral  (and  possibly  ambient)  gas  by  a  translating  FRC 
plasmoid.  The  main  idea  of  such  an  entrainment  is  the  interaction  of  fast  ions  with  slow  neutrals,  subsequent 
charge  exchange  reactions,  and  the  creation  of  slow  ions  and  fast  neutrals.  Fast  neutrals  will  then  exit  the 
thruster  and  produce  thrust,  while  slow  ions  may  again  be  accelerated.  The  repetition  of  this  process 
will  allow  drastic  reduction  of  energy  loss  on  ionization,  while  gaining  thrust  through  multiple  collisions  of 
accelerated  ions  with  slow  neutrals.  A  significant  problem  that  may  reduce  the  usefulness  of  the  neutral 
entrainment  is  the  electron  impact  ionization  reactions.  These  reactions  may  not  just  reduce  the  benefits  of 
entrainment,  but  potentially  diminish  its  usefulness. 

In  recent  paper,7  chemical  reactions  between  ions  and  electrons  of  a  translated  at  high  velocities  field  re¬ 
versed  plasmoid  and  neutral  gas  particles  were  examined.  Importance  of  different  reaction  mechanisms,  such 
as  charge  exchange  and  electron  impact  ionization,  was  analyzed  for  conditions  typical  for  FRC  thrusters. 
Then,  an  implicit  Particle-In-Cell  code,  Celeste3D,8  extended  to  include  plasma-neutral  and  neutral-neutral 
elastic  and  inelastic  processes  and  Coulomb  collisions,  was  used  to  model  a  planar  flow  between  an  FRC 
plasmoid  initially  at  a  Schmid-Burgk  equilibrium,  and  neutral  gas  with  varying  relative  velocities,  plasma 
temperatures,  and  neutral  gas  densities.  The  contribution  of  charge  exchange  and  electron  impact  ionization 
reactions  was  analyzed,  and  the  sensitivity  of  the  entrainment  efficiency  to  the  plasmoid  translation  velocity 
and  neutral  density  was  evaluated. 

In  the  present  work,  which  is  the  continuation  of  our  previous  effort,7  the  focus  is  on  (i)  implementing  the 
electronic  excitation  processes  and  (ii)  analysis  of  thruster  efficiencies  for  a  neon  FRC  plasmoid  with  various 
plasma  temperatures.  Efficiency  of  using  local  gas  injection  is  evaluated,  and  the  limits  of  thrust  increase 
through  the  neutral  density  variation  are  examined.  In  the  next  section,  the  proposed  model  for  collisional 
radiative  processes  in  neutral  entrainment  of  a  neon  plasmoid  is  presented.  Then,  the  collisional  relaxation 
in  a  spatially  homogeneous  thermal  bath  is  studied.  Finally,  the  implementation  of  collision  processes  to 
Celeste3D  and  modeling  of  a  2D  FRC  neutral  entrainment  is  discussed. 

II.  Collisional  radiative  and  reacting  processes  in  neon  neutral  entrainment 

The  main  focus  of  this  work  is  on  neon  propellant,  since  it  appears  to  benefit  from  high  ionization 
threshold,  and  thus  relatively  low  losses  due  to  ionization  in  collisions  between  hot  electrons  of  the  plasmoid 
and  cold  neutrals  of  the  entraining  gas.7  Xenon,  which  may  be  attractive  due  to  its  high  mass,  has  much 
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lower  ionization  threshold,  and  the  associated  energy  losses  may  significantly  decrease  the  entrainment 
efficiency.  Helium,  conversely,  is  difficult  to  ionize  but  its  low  mass  make  it  poor  propellant.  In  the  first 
work  where  kinetic  modeling  was  used  to  analyze  the  plasmoid  translation  and  neutral  entrainment,7  the 
reaction  processes  between  a  neon  plasmoid  and  neutral  gas  included  only  the  charge  exchange  and  electron 
impact  ionization  reaction.  The  electronic  excitation,  and  the  resulting  loss  of  plasmoid  energy  to  radiation, 
was  not  considered,  even  though  it  is  an  important  process  with  a  rate  exceeding  that  of  the  electron  impact 
ionization  at  plasma  temperatures  of  about  or  below  5  eV.  In  the  present  work,  all  key  relaxation  and  reaction 
processes  are  included  in  the  model. 

The  following  neon  states  are  known  to  exist  in  a  neon  plasma: 

1.  Neutral  atoms  in  the  ground  electronic  state, 

2.  Neutral  atoms  in  nretastable  states  with  lifetimes  in  tens  of  seconds  (effectively  infinite  on  a  sub- 
millisecond  scale  of  the  FRC  entrainment  process), 

3.  Neutral  atoms  in  resonant  states  that  quickly  decay  into  ground  state  (such  a  decay  is  typically  fast 
enough  to  be  considered  instantaneous), 

4.  Neutral  atoms  in  non-resonant  excited  states  with  short  radiative  lifetimes  (14.5-27  ns).9  From  these 
states,  the  atoms  may  radiatively  decay  either  to  a  metastable  or  a  resonant  state. 

5.  Positively  charged  ions  (only  singly  charged  ions  need  to  be  considered  in  a  low -Z  FRC  thruster 
environment). 

The  lifetimes  of  resonant  and  non-resonant  states  are  typically  in  tens  of  nanoseconds,  while  the  charac¬ 
teristic  time  of  neutral  entrainment,  which  is  the  time  it  takes  an  FRC  plasmoid  moving  at  about  20  knr/s 
to  pass  through  the  neutral  field,  is  in  tens  of  microseconds.  Therefore,  it  may  be  assumed  without  any  loss 
of  accuracy  that  the  resonant  and  non-resonant  states  emit  photons  immediately  after  the  corresponding 
collision,  and  thus  the  collision/reaction  process  is  reduced  to  the  energy  loss  of  the  colliding  particle  pair. 
That  means  that  in  the  kinetic  model,  only  three  states  need  to  be  considered:  neutral  atoms,  metastable 
atoms,  and  ions.  Since  neon  has  two  nretastable  levels  with  closely  spaced  energies  (16.619  eV  and  16.716 
eV),  those  can  be  effectively  lumped  together  without  any  noticeable  loss  of  accuracy.  The  resulting  model 
appears  as  follows. 

Part  I:  Elastic  collision  processes 

•  Neon  ground  state  atom  elastic  collisions  with  electrons.  The  cross-sections  for  high  energy  elastic 
collisions  are  taken  from  Ref.  10,  11,  and  the  values  for  low  impact  energies  are  taken  from  Ref.  12. 

•  Elastic  collisions  between  neutral  atoms.  The  variable  hard  sphere  model  with  parameters  from  Ref.  13 
is  applied  for  this  interaction. 

•  Coulomb  collisions  between  charged  particles.  The  model14  is  implemented  to  compute  these  collisions. 

Part  II:  Plasmachemical  reactions: 

•  Charge  exchange  reactions.  Similar  to  the  previous  work,7  the  charge  exchange  reaction  cross  section 
of  Ref.  15  is  used  here.  The  after-collision  velocities  of  the  colliding  ion  and  neutral  are  swapped. 

•  Elastic  collisions  between  neutrals  and  electrons.  The  cross  section  for  this  process  was  assumed  to  be 
one-half  of  the  charge  exchange  cross  section,  according  to  the  recommendation  of  Ref.  15. 

Part  III:  Collision  radiative  processes 

•  Electron  impact  ionization  of  ground  state  atoms,  Ne  +e  — >  Ne+  +  e  +  e.  For  this  reaction,  the  cross 
sections  from  SIGLO  database16  are  used  in  the  computations. 

•  Electron  impact  excitation  of  ground  state  atoms  to  metastable  states,  Ne+e  — ►  Ne*M  +  e.  The  cross- 
sections  for  this  process  are  taken  from  Ref.  17. 
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•  Electron  impact  excitation  from  the  ground  state  to  a  short-lived  excited  state,  followed  by  radiative 
decay  into  either  metastable  or  resonance  state.  Since  the  resonance  state  is  short-lived,  and  the  cross- 
section  for  the  excitation  to  nretastable  state  already  takes  into  account  the  cascade  contributions  from 
higher  levels,  the  product  of  this  reaction  should  be  a  ground-state  atom,  Ne  +e  — >  Ne+e  +  hv.  In  this 
work,  the  cross  sections  recommended  in  Ref.  18  are  used. 

•  Electron  impact  excitation  of  ground  state  atoms  to  resonant  states,  followed  by  the  radiative  decay 
to  the  ground  state,  Ne  +e  —>  Ne+e  +  hv.  For  this  process,  the  cross-section  are  taken  from  Ref.  19. 

•  Transition  of  metastable  atoms  to  resonant  states,  with  the  subsequent  de-excitation,  Ne^+e  — >Ne+e+ 
hv.  The  rate  constant  for  this  process  is  equal  to  2  •  10“  7  cnr3/s,20  and  the  cross-sections  for  kinetic 
modeling  are  calculated  using  the  total  collision  energy  model.21 

•  Electron  impact  ionization  of  metastable  atoms,  Ne^  +  e  — >-Ne+  +  e  +  e.  The  cross  sections  for  this 
process  are  taken  from  Ref.  22. 

•  Electron  impact  excitation  of  metastable  atoms  to  non-resonant  exited  states,  Ne/f  +  e  — >  Ne*  +  e. 
Since  the  products  of  this  reaction  are  short-lived,  the  reaction  product  is  a  ground-state  atom.  The 
reaction  cross-sections  follow  from  Ref.  18. 

Due  to  the  lack  of  detailed,  differential  collision  cross  sections  for  most  of  the  above  processes,  an  isotropic 
scattering  is  assumed  to  calculate  after-collision  velocities  of  atoms  and  electrons  in  all  interactions  with 
the  exception  of  Coulomb  collisions  and  charge  exchange  reactions.  This  assumption  is  not  expected  to 
significantly  impact  the  results  presented  in  the  following  sections. 

III.  Collisional  relaxation  of  plasma  and  neutral  gas  in  OD  heat  bath 

Consider  first  the  impact  of  the  plasmoid  temperature  on  the  time  evolution  of  plasma  and  neutral 
properties  in  an  adiabatic  heat  bath  with  key  macroscopic  properties,  such  as  relative  velocity  between 
neutrals  and  plasma,  and  plasma  temperature,  that  correspond  to  conditions  typical  for  an  FRC  thruster. 
It  is  assumed  that  the  neutral  and  charged  species,  initially  at  equilibrium  at  their  respective  temperatures, 
have  a  relative  velocity  of  20  km/s  between  them;  due  to  the  symmetry  of  the  problem,  the  average  plasma 
velocity  of  0  and  neutral  velocity  of  20  km/s  is  assumed.  The  initial  neutral  temperature  is  10  K,  and  the 
plasma  temperature  is  either  5  eV  or  10  eV;  both  neutral  and  plasma  density  are  1018  m-3.  To  simulate 
the  neutral  entrainment  process,  a  0D  direct  simulation  Monte  Carlo  (DSMC)  code  was  developed  that 
incorporates  all  processes  listed  in  the  previous  section. 

The  temperature  profiles  of  neutrals  (T„),  ions  (T)),  and  electrons  (Te)  as  a  function  of  time  are  shown 
in  Fig.  2  (left)  for  an  initial  plasma  temperature  of  5  eV.  During  the  first  few  microseconds,  there  are  not 
enough  ionization  reactions  for  the  electron  temperature  to  noticeably  change.  At  the  same  time,  the  high- 
energy  elastic  collisions  between  neutrals  and  ions  (their  relative  collision  velocity  of  20  km/s  corresponds  to 
an  energy  of  20  eV,  and  thus  a  temperature  of  about  10  eV)  result  in  the  increase  in  both  ion  and  neutral 
temperature.  The  maximum  ion  and  neutral  temperature  of  about  9  eV  is  reached  at  about  100  /is.  Note 
that  the  maximum  value  is  impacted  by  the  ionization  and  radiation  losses,  that  start  to  contribute  after 
about  10  /is.  The  electron  temperature  does  not  have  a  visible  maximum  due  to  relatively  slow  relaxation 
on  heavy  particles,  and  also  due  to  energy  losses  to  the  electronic  excitation  /  radiation  and  ionization.  The 
relative  contribution  of  the  ionization  and  the  collisional  radiative  processes  is  illustrated  in  Fig.  2  through  the 
simulation  of  the  plasma/neutral  interaction  with  only  the  charge  exchange  and  electron  impact  ionization 
included  (the  corresponding  electron  temperature  is  labeled  “Te  no  CR”).  It  is  clear  that  the  collisional 
radiation  is  more  important  for  the  relatively  low  temperature  of  5  eV  than  the  ionization  processes,  as  it  is 
responsible  for  over  70%  of  the  electron  temperature  decrease.  By  the  time  of  1  ms,  about  three  percent  of 
electrons  and  neutrals  have  participated  in  ionization. 

The  temperature  evolution  for  a  hotter  plasma  is  qualitatively  similar  to  5  eV,  but  the  impact  of  the 
radiation  and  ionization  is  much  more  pronounced  in  this  case  due  to  the  increased  rates  of  these  processes. 
After  1  ms,  the  electron  temperature  decreases  by  almost  60%,  as  about  15%  of  electrons  and  neutrals  par¬ 
ticipated  in  electron  impact  ionization..  Most  of  that  decrease  is  due  to  the  ionization,  as  the  corresponding 
temperature  for  the  no-excitation  case  is  about  50%  of  the  initial  temperature.  Similar  to  the  5  eV  case, 
the  neutral  and  ion  temperatures  increase  over  the  first  100  /is,  and  then  decrease  due  to  the  excitation 
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and  ionization.  It  is  important  to  note  that  for  both  initial  plasma  temperatures  there  is  a  strong  thermal 
non-equilibrium,  with  the  species  temperatures  significantly  different  due  to  the  elastic  and  inelastic  collision 
processes. 


Figure  2.  Species  temperature  relaxation  for  initial  plasma  temperature  of  5  eV  (left)  and  10  eV  (right). 

The  non-equilibrium  is  manifested  not  only  at  the  macroscopic  level  with  the  species  temperature  dif¬ 
ference,  but  also  at  the  microscopic  level  of  the  velocity  distribution  functions.  The  distribution  functions 
of  ion  velocities  in  the  direction  of  neutral-plasma  relative  velocity  at  different  time  moments  are  shown  in 
Fig.  3  (left).  Initially  Maxwellian  (Gaussian)  distribution  becomes  clearly  bimodal  after  the  first  few  mi¬ 
croseconds.  The  peak  at  the  20  km/s  corresponds  to  ions  that  participated  in  charge  exchange  reactions  and 
thus  acquired  the  velocity  of  neutrals.  The  bimodality  essentially  disappears  after  about  50  /is,  but  the  ion 
distribution  function  is  still  visibly  non-Maxwellian  even  after  1  ms.  The  neutral  atom  distribution  function 
(not  shown  here)  is  also  non-equilibrium,  and  qualitatively  similar  to  that  of  ions.  The  non-equilibrium  in  the 
velocity  distribution  functions  shows  the  importance  of  using  a  kinetic  approach  for  modeling  FRC  neutral 
entrainment,  since  a  continuum  approximation  would  typically  undepredict  (and  in  some  cases  over-predict) 
the  actual  reaction  rates. 

The  electron  velocity  distribution  function  is  fairly  close  to  equilibrium  throughout  the  relaxation  process, 
as  illustrated  in  Fig.  3  (right)  for  a  time  moment  of  1  ms.  The  distribution  function  profile  is  nearly 
Maxwellian  both  for  the  baseline  model  that  includes  all  interaction  types,  and  the  model  that  does  not 
take  into  account  the  excitation  processes.  The  change  in  shape  between  the  two  models  corresponds  to 
the  temperature  difference  between  them  discussed  earlier  (see  fig.  2,  left).  It  is  important  to  note  that 
the  principal  mechanism  responsible  for  the  equilibration  of  the  electron  velocity  distribution  function  is 
electron-electron  Coulomb  collisions.  The  distribution  function  obtained  in  a  simulation  that  does  not 
include  these  collisions  is  also  presented  in  Fig.  3  (right);  it  shows  a  significant  depletion  of  high  velocity  tails 
of  the  distribution  function.  As  a  result,  ignoring  Coulomb  collisions  would  strongly  decrease  the  number  of 
electronic  excitations  and  electron  impact  ionizations. 


IV.  Numerical  approach  and  flow  conditions  for  2D  FRC  entrainment 

Modeling  of  the  interaction  between  an  FRC  plasmoid  and  the  entraining  gas  is  conducted  in  this  work 
with  an  implicit  three-dimensional  PIC  code  Celeste3D  that  solves  the  full  set  of  Maxwell- Vlasov  equation 
and  has  been  extensively  used  in  the  past  for  various  astropliysical  and  laboratory  plasma  problems. 23:24  The 
implicit  moment  formulation  of  the  PIC  method  implemented  in  Celeste  results  in  highly  efficient  simulations 
based  on  ion  length  and  time  scales,  and  not  electron  scales  as  explicit  methods  do,  while  retaining  the  kinetic 
effects  of  both  the  electrons  and  the  ions.  An  explicit  simulation  requires  the  time  step  to  be  At  <  2/ojpe, 
and  the  spatial  cell  size  to  be  Ax  <  £Ae  in  order  to  avoid  the  finite  grid  instability.  Here,  uipe  is  the  electron 
plasma  frequency,  and  Xe  is  the  electron  Debye  length.  In  an  implicit  simulation,  these  requirements  are 
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Figure  3.  Ion  velocity  distribution  function  as  a  function  of  time  (left).  Electron  velocity  distribution  function 
at  1  ms  (right).  The  initial  plasma  temperature  of  5  eV. 


replaced  by  an  accuracy  condition  related  to  the  conservation  of  energy,  At  <  fracAxce,  where  ce  is  the 
electron  thermal  speed. 

The  original  Celeste3D  models  the  evolution  of  ion  and  electron  populations  as  described  by  the  Vlasov 
equation,  coupled  with  the  solution  of  full  Maxwell  equations;  it  does  not  simulate  neutral  particles.  To  study 
the  neutral  entrainment  process,  Celeste3D  was  extended  in  the  previous  work'  to  include  neutral  transport 
and  collisional  relaxation.  A  Direct  Simulation  Monte  Carlo  based  capability  has  been  added  to  Celeste  that 
models  neutral  transport  and  collisions.  The  collision  processes  included  in  the  present  model  are  discussed 
in  Section  II.  As  all  species  have  different  weights,  a  weighting  scheme25  is  applied.  The  majorant  collision 
frequency26  of  the  DSMC  method  is  used  to  simulate  the  collision  process  in  cells. 

Although  the  actual  entrainment  geometry  is  three-dimensional,  as  a  nearly  axisymmetric  plasmoid 
interacts  with  neutral  gas  supplied  from  two  or  more  azimuthal  injectors,  three-dimensional  PIC  modeling 
is  extremely  time  consuming,  and  in  this  work  it  was  replaced  by  a  2D  (planar  flow).  The  result  is  a  much 
larger  number  of  particles  per  cell  and  thus  smaller  statistical  scatter.  The  grid  was  40x80,  and  the  number 
of  simulated  ions,  electrons,  and  neutrals  per  cell  was  approximately  64,  100,  and  512,  respectively.  Since  the 
cost  of  an  implicit  simulation  is  a  direct  function  of  the  electron  mass  (see  above,  the  computational  time  is 
inversely  proportional  to  the  square  root  of  the  electron  mass  me),  the  simulation  efficiency  may  be  further 
enhanced  through  the  introduction  of  the  weighted  electron  mass,  m'e  =  Wme,  where  the  constant  W  >  1. 
The  ion-to-electron  mass  ratio  has  to  be  high  enough  to  preserve  the  kinetic  effects,  and  the  value  of  rrii/m'e 
on  the  order  of  100  is  usually  sufficient.23  In  this  work,  m,i/m'e  =  100  is  used  (modeling  of  the  baseline  case 
of  a  5  eV  plasmoid  with  a  neutral  entrainment  for  a  ratio  of  rrii/m'e  =  300  have  shown  identical  results).  The 
initial  conditions  used  to  set  an  equilibrium  plasmoid  is  the  a  Schmid-Burgk  equilibrium27  with  the  average 
plasma  density  of  1018  m-3  and  two  plasma  temperatures  of  5  and  10  eV.  In  all  computations,  the  size  of 
the  computational  domain 

To  simplify  the  inflow  boundary  conditions,  the  neutral-plasma  interaction  is  examined  in  the  reference 
frame  of  the  plasmoid,  so  that  the  neutral  gas  is  injected  into  the  computational  domain  from  the  left 
boundary,  and  then  passes  through  the  domain  from  left  to  right.  The  properties  of  the  injected  neutral 
gas  are  changed  to  study  their  impact  on  the  neutral-plasma  interaction.  The  difference  from  the  previous 
work,7  in  addition  to  a  much  more  complete  collisional  relaxation  model,  is  the  local  injection  of  neutrals. 
While  in  the  previous  work  the  neutrals  were  introduced  from  the  entire  left  boundary,  here  the  neutrals  are 
injected  only  in  the  central  part  of  the  flow  (the  injection  area  was  0.02  m).  This  allows  one  to  to  study  the 
trajectories  of  the  after-collision  atoms,  and  analyze  the  usability  of  local  injection  concept. 

In  all  computations,  open  boundary  conditions  for  neutrals  and  plasma  were  imposed  at  the  left  and 
right  boundaries,  and  a  conducting  wall  with  specular  reflection  was  set  at  the  top  and  bottom  boundaries. 
Note  that  computations  were  also  performed  with  no  neutral  inflow  to  analyze  the  temporal  evolution  of 
the  plasmoid;  it  was  observed  that  the  equilibrium  is  maintained  for  at  least  10,000  ion  plasma  periods. 
This  amounts  to  over  30  [is  for  the  conditions  under  consideration,  and  is  longer  than  the  plasmoid-neutral 
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interaction  time.  The  macroparameters  presented  below  are  the  result  of  an  instantaneous  average  (ensemble 
sampling)  of  gas  and  plasma  properties  over  cells.  The  distribution  functions  of  neutral  atom  velocities  were 
also  ensemble  averaged.  The  thrust  force  was  calculated  in  this  work  by  an  integration  over  a  cross  section 
that  moves  with  a  velocity  equal  to  that  of  the  neutral  gas.  An  axial  symmetry  around  the  lower  boundary 
is  assumed  in  the  integration  in  order  to  provide  the  link  between  the  computed  force  and  the  actual  FRC 
thruster. 

The  schematics  of  the  numerical  setup  is  illustrated  in  Fig.  4  where  the  ion  and  neutral  densities  are 
shown  at  time  t  =  8  /zs.  Hereafter,  t  =  0  corresponds  to  the  initial  state  when  the  neutral  gas  starts  to  flow 
into  the  computational  domain.  At  t  =  8  /us,  the  injected  at  X  =  0  within  0.09  <  Y  <  0.11  m  neutrals  have 
crossed  slightly  more  than  half  of  the  domain.  The  neutral  bulk  velocity  in  X  direction  is  20  krn/s,  much 
higher  than  the  average  thermal  velocity  of  about  0.1  krn/s,  and  thus  the  divergence  of  the  flow  is  fairly 
small.  Neutrals  that  appear  outside  the  0.09  <  Y  <  0.11  m  region  have  traveled  there  after  being  scattered 
as  a  result  of  ion-neutral  collisions.  The  ion  density  is  maximum  at  the  center  of  the  domain,  and  that  is 
the  region  that  affects  neutrals  most. 


X,  m 


Figure  4.  FRC  neutral  entrainment  setup.  Field:  ion  number  density  (10ls  m  3),  isolines:  neutral  gas  number 
density  (1018  m~3). 


V.  FRC  neutral  entrainment  results 

The  time-dependent  evolution  of  neutral  gas  penetrating  the  plasmoid  is  given  in  Fig.  5,  left,  where  the 
gas  density  is  shown  for  two  initial  plasma  temperatures  of  5  and  10  eV  and  three  time  moments,  t  =  0.08,  8, 
and  16  /is.  After  the  first  microsecond,  there  is  no  visible  change  of  neutral  gas  density  from  its  free  stream 
value.  As  neutral  gas  starts  to  penetrate  denser  plasma  region,  some  scattering  of  neutrals  to  the  side  is 
observed.  However,  after  the  neutral  jet  passed  through  the  plasmoid  center,  there  is  still  no  neutrals  at 
the  walls  of  the  FRC  chamber,  primarily  because  there  was  not  enough  time  for  reflected  neutrals  to  reach 
them.  The  average  thermal  velocity  of  neutrals  after  collisions  with  ions  corresponds  to  the  temperature 
of  the  plasmoid,  and  this  is  about  7  km/s  for  5  eV  and  about  10  km/s  for  10  eV.  After  the  neutral  gas 
passed  through  the  plasmoid  ( t  =  16  /zs),  the  reflected  neutrals  reach  the  chamber  walls,  but  there  density 
in  that  region  is  fairly  small,  up  to  three  orders  of  magnitude  smaller  than  that  in  the  coreflow.  This  is 
a  combined  effect  of  relatively  small  number  of  collisions  that  mostly  occur  in  the  center  of  the  plasmoid, 
lower  after-collision  neutral  velocities,  and  the  preferential  direction  of  these  velocities  (aligned  with  neutral 
jet).  Note  that  the  impact  of  plasmoid  temperature  on  neutral  density  is  relatively  small,  as  only  a  few 
percent  smaller  number  of  neutrals  reaches  the  right  boundary  for  10  eV.  This  few  percent  loss  is  primarily 
attributed  to  electron  impact  ionization  reactions. 

The  evolution  of  ion  number  density  for  two  plasmoid  temperatures  as  the  result  of  the  neutral-plasma 
interaction  is  shown  in  Fig.  5,  right.  The  change  in  ion  density  during  the  first  8  /zs  is  small,  and  related 
mostly  to  the  statistical  perturbations  and  some  minor  change  in  density  distribution  as  the  plasmoid  slightly 
deviates  from  the  Schmid-Burgk  equilibrium,  and  not  to  the  entrainment  process.  The  only  neutral-related 
change  at  that  time  is  noticeable  in  the  plasmoid  coreflow.  The  change  between  the  8  /zs  and  16  /zs  fields 
is  more  significant.  As  expected,  it  is  observed  only  in  the  coreflow,  the  region  where  the  neutral  jet  passes 
through.  At  5  eV,  the  number  of  ionization  reaction  is  much  smaller  than  at  10  eV,  and  thus  the  neutral 
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entrainment  is  less  pronounced.  At  10  eV,  the  ion  density  in  the  coreflow  increases  by  up  to  10%  with  respect 
to  the  initial  equilibrium  values.  There  is  also  some  momentum  transfer  from  neutral  jet  to  the  plasmoid, 
as  the  latter  visibly  shifts  to  the  right  at  t  =  16  /xs. 
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Figure  5.  Neon  number  density  (m  3),  left,  and  ion  density  (m  3), right,  for  5  eV  (left  columns)  and  10  eV 
(right  columns)  at  different  time  moments  since  the  injection.  Injected  neutral  density  is  1018  atom/ml 


As  follows  from  Fig.  5,  left,  the  number  of  neutrals  that  reach  the  right  boundary  of  the  computational 
domain  for  the  neutral  jet  density  of  1018  atom/m3  is  about  87%  and  85%  for  5  eV  and  10  eV  plasmoids, 
respectively.  Since  the  number  density  of  neutrals  is  close  to  the  average  number  density  of  the  plasmoid, 
that  implies  that  only  on  the  order  of  10%  of  the  plasma  particles  had  collisions  with  neutrals.  There  are  two 
obvious  ways  to  increase  that  number:  (i)  increase  the  interaction  time,  or  (ii)  increase  neutral  density.  The 
first  approach  is  not  efficient  since  it  requires  either  decreasing  the  plasmoid  velocity,  and  thus  will  reduce 
the  thruster  specific  impulse,  or  increasing  the  thruster  size.  The  second  approach  may  be  more  beneficial, 
especially  remembering  that  in  the  actual  thruster,  the  pulsed  mode  of  the  fast  plasmoid  along  with  low 
velocities  of  neutrals  may  allow  for  a  multi-pass  entrainment  of  neutrals.  In  order  to  analyze  the  impact  of 
higher  neutral  density,  the  computations  were  conducted  for  densities  that  are  5  and  25  times  higher  than 
the  baseline  of  n  =  1018  atom/m3. 

The  neutral  density  field  for  n  =  2.5  x  1019  atom/m3  at  different  time  moments  is  shown  in  Fig.  6,  left, 
for  the  two  plasmoid  temperatures  under  consideration.  Qualitatively,  the  plot  looks  similar  to  that  for  the 
25  time  lower  neutral  density.  No  significant  deflection  is  noticeable  after  the  first  microsecond  after  the 
injection;  at  8  /x s,  there  is  a  noticeable  amount  of  neutrals  scatter  to  the  side,  although  none  of  those  have 
reached  the  chamber  walls  yet.  At  16  /us,  the  neutrals  reflected  after  collisions  with  ions  reach  the  chamber, 
although  their  density  near  the  wall  is  full  three  orders  of  magnitude  smaller  than  in  the  coreflow.  More 
importantly,  almost  12%  of  neutral  particles  collide  with  plasma  particles  as  the  neutral  jet  crosses  the  5  eV 
plasmoid;  for  a  10  eV  plasmoid,  this  number  reaches  15%,  similar  to  the  lower  density  case. 

Since  the  neutral  density  is  much  higher  than  that  of  the  plasmoid,  one  can  expect  significant  entrainment 
in  this  case.  And  this  is  indeed  the  case,  as  Fig.  6,  right,  clearly  shows.  There  is  a  significant  increase  of 
plasma  density  even  at  t  =  8  /is.  At  that  time,  the  ion  and  electron  density  in  the  center  of  the  plasmoid  is 
over  50%  higher  than  the  corresponding  equilibrium  value  for  the  5  eV  plasmoid,  and  is  almost  two  times 
higher  for  the  higher  temperature  of  10  eV.  Note  that  such  a  difference  between  5  eV  and  10  eV  is  consistent 
with  the  temperature  dependence  of  the  electron  impact  ionization  reaction. '  Further  interaction  between 
the  plasmoid  and  neutral  gas  ( t  =  16  /is)  results  in  a  lower  maximum  plasma  density  that  for  t  =  8  /is. 
This  is  related  primarily  to  the  neutral-to-plasma  momentum  transfer.  At  t  =  16  fi s,  the  plasmoid  loses  its 
equilibrium  ellipsoidal  shape,  and  significant  part  of  it  translates  with  the  neutrals.  This  is  due  to  charge 
exchange  and,  to  a  lesser  extent,  elastic  collisions,  both  of  with  transfer  the  X  momentum  to  ions.  While 
the  total  number  of  plasma  particles  is  significantly  larger  at  16  /is  than  at  8  /xs,  this  number  is  spread  over 
significant  distance. 

The  neutral  entrainment  also  results  in  significant  increase  of  ion  temperature,  discussed  in  detail  in  the 
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Figure  6.  Same  as  Fig.  5,  for  neutral  density  of  25  x  1018  atom/m'1. 


previous  section.  Due  to  highly  energetic  elastic  collisions  between  neutrals  and  ions  (the  average  collision 
energy  is  about  40  eV),  as  well  as  charge  exchange  reactions  that  increase  ion  velocities  to  20  krn/s,  the  ion 
temperature  increases  to  almost  10  eV  when  the  initial  plasmoid  temperature  was  5  eV,  and  to  over  13  eV  for 
the  10  eV  plasmoid,  as  shown  in  Fig.  7,  left.  This  result  is  consistent  with  the  earlier  Fig.  2.  Note  also  that 
at  16  /x s  the  high  ion  temperature  region  splits  in  two  parts,  above  and  below  the  domain  symmetry  axis, 
which  is  attributed  to  the  plasmoid  translation  and  ion  mixing  in  these  regions.  The  computations  show  that 
the  electron  temperature  (not  shown  here)  is  practically  not  affected  by  the  plasmoid  entrainment  during 
the  plasma-neutral  interaction.  Due  to  the  low  rate  of  momentum  and  energy  transfer  in  the  electron-heavy 
particle  collisions,  the  electron  temperature  does  not  change  in  the  computations  by  more  than  a  few  percent 
(decreases  by  about  5%  for  a  10  eV  plasmoid).  This  result  is  also  in  line  with  what  was  observed  earlier  in 
the  homogeneous  heat  bath  relaxation. 

All  plasma  and  neutral  species  are  in  thermal  nonequilibrium,  which  is  manifested  in  differences  in  their 
translational  temperatures  as  well  as  in  non-Maxwellian  distribution  functions.  The  elector  and  ion  velocity 
distribution  in  the  2D  FRC  entrainment  are  similar  to  those  considered  in  the  previous  section,  so  here  we 
only  present  the  distribution  function  of  X  component  velocity  for  neutral  atoms.  The  velocity  distribution 
function  of  neutrals  is  shown  in  Fig.  7,  right.  Since  single-cell  transient  results  have  considerable  noise,  the 
average  of  all  cells  is  presented  in  this  figure.  Note  that  the  velocity  distributions  weakly  depend  on  the 
neutral  gas  density,  and  thus  only  the  case  of  1018  atom/m3  is  given  here.  As  clearly  seen  in  Fig.  7,  right,  the 
neutral  velocity  is  bimodal,  with  a  sharp  peak  at  -20  km/s  corresponding  to  the  injected  atoms  that  did  not 
have  any  collisions  with  plasma,  and  the  second,  much  lower  maximum,  at  U  =  0  corresponding  to  the  atoms 
collided  with  ions.  Fairly  quick  equilibration  after  collisions  with  the  isotopic  scattering  results  in  the  second 
maximum  having  a  shape  of  a  Maxwellian  distribution  with  the  corresponding  plasma  temperature.  This  is 
illustrated  by  comparison  with  the  Mott-Smith  bimodal  solution  with  the  mode  temperatures  of  10  K  and 
lOeV.  The  ratio  of  collided  to  non-collided  neutrals  is  approximately  0.1  for  both  plasmoid  temperatures. 

Consider  now  the  contribution  of  different  gas  and  plasma  species  to  the  thrust  force,  illustrated  in 
Fig.  8  for  a  5  eV  plasmoid  and  a  neutral  density  of  1018  atom/m3.  In  this  figure,  p  denotes  the  contribution 
of  the  pressure  component  (thermal  motion  of  particles)  of  the  corresponding  species,  and  pU2  denotes 
the  contribution  of  the  momentum  component  (translational  motion).  The  momentum  component  of  the 
electrons,  as  well  as  pressure  component  of  the  neutral  gas,  are  negligible  and  thus  not  plotted  here.  The 
next  important  contributors  are  ion  and  electron  pressures,  which  almost  coincide  in  the  considered  case; 
their  combined  effect  amounts  at  its  peak  to  about  10%  of  the  total  thrust  force. 

It  is  important  to  remember  that  there  are  several  key  differences  of  the  considered  setup  from  an  actual 
operating  FRC  thruster.  First,  the  reference  frame  of  the  plasmoid  is  chosen  in  the  computations.  Second, 
the  numerical  simulation  is  two-dimensional.  Third,  and  probably  most  important,  the  modeling  lacks  an 
accelerating  magnetic  field.  Such  a  field  would  transfer  momentum  to  the  plasmoid,  and  replenish  the  mo¬ 
mentum  loss  caused  by  the  collisions  with  the  neutrals.  As  the  result  of  the  latter  numerical  approximation, 
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Figure  7.  Left:  ion  translational  temperature  at  different  time  moments  for  neutral  density  of  25xl018  atom/m3. 
Right:  neutral  velocity  distribution  function  at  t  =  16  ps  for  neutral  density  of  1018  atom/m3. 


and  due  to  the  conservation  of  mass,  momentum,  and  energy  in  the  numerical  modeling,  the  total  thrust 
force  computed  over  neutral  and  plasma  species  does  not  change  when  the  physical  model  or  flow  conditions 
are  modified.  The  indication  of  the  entrainment  efficiency  is  therefore  the  amount  of  energy  and  momentum 
transferred  from  plasma  (primarily  ion)  species  to  the  neutrals.  The  larger  the  decrease  in  thrust  force 
produced  by  the  ions,  the  stronger  the  entrainment.  The  situation  is  somewhat  complicated  by  the  presence 
of  several  mechanisms  that  impact  the  ion  component  of  thrust,  primarily  the  charge  exchange  and  elastic 
ion-neutral  collisions,  and  electron  impact  ionization  of  neutrals,  but  nonetheless  the  decrease  in  the  ion 
thrust  component  is  the  key  factor  of  the  entrainment,  and  thus  considered  below. 

Figure  8,  right,  presents  the  temporal  change  in  thrust  for  three  values  of  neutral  density  considered  in 
this  work.  The  results  show  that  increase  in  neutral  density  leads  to  considerable  decrease  in  ion  contribution 
to  thrust,  which  decreases  at  its  peak  by  about  10%  as  the  neutral  density  changes  from  1018  atom/m3  to 
5  x  1018  atom/m3,  and  another  25%  as  density  goes  up  to  25  x  1018  atom/m3.  The  decrease  is  less  then 
linear  due  to  a  non-linear  momentum  and  mass  transfer  during  the  entrainment,  the  threshold  nature  of  the 
ionization  and  collisional  radiative  processes,  and  the  non-linear  charge  exchange  rate  that  is  impacted  by 
temperature  variations  shown  in  Fig.  7,  left.  Note  that  the  introduction  of  neutrals  only  in  the  coreflow, 
as  compared  to  the  uniform  injection  across  the  chamber,  is  beneficial  only  as  long  as  the  plasma  density 
is  higher  in  the  center,  since  the  entrainment  efficiency  is  governed  both  by  neutral  and  plasma  density. 
The  results  for  a  10  eV  plasma  (not  shown  here)  are  quantitatively  similar  to  the  5  eV  case  for  all  neutral 
densities  under  consideration. 


VI.  Conclusions 

The  paper  presents  the  results  of  numerical  modeling  of  neutral  entrainment  in  a  Field  Reversed  Config¬ 
uration  thruster.  The  numerical  approach  chosen  is  the  combination  of  the  implicit  Particle-In-Cell  method 
for  modeling  plasma  and  the  direct  simulation  Monte  Carlo  method  for  modeling  neutral  particles.  The 
multi-dimensional  PIC  code  Celeste3D  extended  with  a  DSMC  capability  was  used  in  this  work  to  model 
a  neon-based  FRC  plasmoid  interacting  with  neon  neutral  gas.  The  modeling  capability  was  extended  as 
compared  to  the  previous  work  to  include  all  relevant  collisional  radiative  processes  in  neon.  The  plasma 
density  of  1018  m-3  and  temperatures  of  5  and  10  eV  were  assumed,  with  the  neutral  density  varied  between 
1018  m-3  and  2.5  x  1019  m~3.  The  computations  were  conducted  in  the  reference  frame  of  the  plasmoid, 
with  neutral  gas  injected  in  the  plasmoid  coreflow  with  a  relative  velocity  of  20  km/s.  The  results  cover 
macroparameters  and  distribution  function,  as  well  as  thruster  performance  related  properties. 

The  numerical  modeling  included  spatially  uniform  thermochemical  relaxation  (adiabatic  heat  bath) 
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Figure  8. 

with  neutral  and  plasma  properties  close  to  those  of  a  typical  FRC  thruster.  The  heat  bath  computations 
showed  that  the  relaxation  process  proceeds  under  conditions  of  strong  thermal  and  chemical  non-equilibrium; 
ion,  electron,  and  neutral  temperatures  strongly  differ,  and  the  ion  and  neutral  gas  distribution  function 
is  strongly  non-Maxwellian.  The  electron  velocity  distribution  function  is  close  to  Maxwellian  when  the 
Coulomb  collisions  are  taken  into  consideration;  otherwise,  electron  impact  ionization  was  found  to  deplete 
high  velocity  tail  of  the  electron  distribution  function.  Electronic  excitation  and  radiation  processes  were 
found  to  noticeably  decrease  primarily  electron  temperature,  and  the  change  is  especially  noticeable  for 
interaction  times  exceeding  10/is. 

Two-dimensional  modeling  of  the  FRC  neutral  entrainment  showed  strong  entrainment  for  neutral  den¬ 
sities  exceeding  5  x  10ls  m-3,  with  momentum  transfer  and  thruster  efficiency  being  non-linear  functions  of 
neutral  density.  The  probability  of  neutral  atom  collision  with  the  plasmoid  was  estimated  at  approximately 
10%,  with  most  of  the  collisions  being  charge  exchange  events.  Electron  impact  ionization  is  not  significant 
for  a  5  eV  plasmoid,  but  becomes  a  factor  at  10  eV.  For  neutral  gas  injected  in  a  10%  area  of  the  chamber, 
the  ion  contribution  to  thrust  was  found  to  decrease  by  up  to  30%  due  to  the  plasma-neutral  momentum 
transfer  when  neutral  density  was  increased  to  2.5  x  1019  m-3. 
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